#include "mesh3d.h"

Mesh3D::Mesh3D(GeoGrid &grid) {
    // 从 GeoGrid 中获取网格信息
    is_expand = grid.is_expand;
    if (is_expand) {
        points = grid.mesh_points_expand;
        tetrahedron_cells = grid.tetrahedron_cells_expand;

        basic_point_size = grid.mesh_points.size();
        basic_tetra_size = grid.tetrahedron_cells.size() / 4;
        expand_point_size = grid.mesh_points_expand.size();
        expand_tetra_size = grid.tetrahedron_cells_expand.size() / 4;
    } else {
        points = grid.mesh_points;
        tetrahedron_cells = grid.tetrahedron_cells;

        basic_point_size = grid.mesh_points.size();
        basic_tetra_size = grid.tetrahedron_cells.size() / 4;
        expand_point_size = basic_point_size;
        expand_tetra_size = basic_tetra_size;
    }

    // initialize data
    young.resize(expand_tetra_size, 0.0f);
    poisson.resize(expand_tetra_size, 0.0f);
    rockden.resize(expand_tetra_size, 0.0f);
    poro.resize(expand_tetra_size, 0.0f);
    if (is_expand) {
        cell_types.resize(expand_tetra_size, ExpandCellType::OTHER); // only for expand mesh
    }
    else {
        cell_types.resize(basic_tetra_size, ExpandCellType::RESERVOIR); // only for basic mesh
    }

    isInterface = grid.mesh_points_type; // only for basic mesh

}

bool Mesh3D::WriteVTK(const std::string filename, bool expand, bool cellData) {
    // 使用 fmt::memory_buffer 将所有数据格式化到内存中
    fmt::memory_buffer buffer;

    // 如果知道大概大小，可以提前 reserve，避免多次分配（可选）
    // buffer.reserve(estimated_size);

    // 写入头部信息
    fmt::format_to(std::back_inserter(buffer),
                   "# vtk DataFile Version 3.0\n"
                   "Mesh\n"
                   "ASCII\n"
                   "DATASET UNSTRUCTURED_GRID\n");

    // 写入点信息
    size_t numPoints = expand ? expand_point_size : basic_point_size;
    fmt::format_to(std::back_inserter(buffer), "POINTS {} float\n", numPoints);
    for (int i = 0; i < numPoints; ++i) {
        fmt::format_to(std::back_inserter(buffer),
                       "{:.10f} {:.10f} {:.10f}\n",
                       points[i].x, points[i].y, points[i].z);
    }

    // 写入单元（cells）
    size_t numCells = expand ? expand_tetra_size : basic_tetra_size;
    // 每个单元前有一个数字 4，加上所有单元的节点索引，总数为 tetrahedron_cells.size() + numCells
    size_t totalCellInts = numCells * (4 + 1);
    fmt::format_to(std::back_inserter(buffer), "CELLS {} {}\n", numCells, totalCellInts);
    for (size_t i = 0; i < 4 * numCells; i += 4) {
        fmt::format_to(std::back_inserter(buffer),
                       "4 {} {} {} {}\n",
                       tetrahedron_cells[i],
                       tetrahedron_cells[i + 1],
                       tetrahedron_cells[i + 2],
                       tetrahedron_cells[i + 3]);
    }

    // 写入单元类型
    fmt::format_to(std::back_inserter(buffer), "CELL_TYPES {}\n", numCells);
    for (size_t i = 0; i < numCells; ++i) {
        fmt::format_to(std::back_inserter(buffer), "10\n");
    }

    // 根据需要写入单元数据（例如 Young 模量、泊松比、岩石密度、孔隙率）
    if (cellData) {
        fmt::format_to(std::back_inserter(buffer), "CELL_DATA {}\n", numCells);

        // Young 模量
        fmt::format_to(std::back_inserter(buffer),
                       "SCALARS young float 1\nLOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i) {
            fmt::format_to(std::back_inserter(buffer), "{:.10f}\n", young[i]);
        }

        // 泊松比
        fmt::format_to(std::back_inserter(buffer),
                       "SCALARS poisson float 1\nLOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i) {
            fmt::format_to(std::back_inserter(buffer), "{:.10f}\n", poisson[i]);
        }

        // 岩石密度
        fmt::format_to(std::back_inserter(buffer),
                       "SCALARS rockden float 1\nLOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i) {
            fmt::format_to(std::back_inserter(buffer), "{:.10f}\n", rockden[i]);
        }

        // 孔隙率
        fmt::format_to(std::back_inserter(buffer),
                       "SCALARS poro float 1\nLOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i) {
            fmt::format_to(std::back_inserter(buffer), "{:.10f}\n", poro[i]);
        }

        // 如果是延拓网格，还需要写入单元类型
        if (is_expand) {
            fmt::format_to(std::back_inserter(buffer),
                           "SCALARS cellType int 1\nLOOKUP_TABLE default\n");
            for (int i = 0; i < numCells; ++i) {
                fmt::format_to(std::back_inserter(buffer), "{}\n", static_cast<int>(cell_types[i]));
            }
        }
    }

    // 写入point data
    if (!isInterface.empty()) {
        fmt::format_to(std::back_inserter(buffer), "POINT_DATA {}\n", numPoints);
        fmt::format_to(std::back_inserter(buffer),
                       "SCALARS isInterface int 1\nLOOKUP_TABLE default\n");
        for (int i = 0; i < numPoints; ++i) {
            if (i < basic_point_size) {
                // if (isInterface[i] == PointType::PT_FAULT || isInterface[i] == PointType::PT_LAYER) {
                //     fmt::format_to(std::back_inserter(buffer), "1\n");
                // } else {
                //     fmt::format_to(std::back_inserter(buffer), "0\n");
                // }
                fmt::format_to(std::back_inserter(buffer), "{}\n", isInterface[i]);
            } else {
                fmt::format_to(std::back_inserter(buffer), "0\n");
            }
        }
    }

    // 将缓冲区内容一次性写入文件
    FILE* vtkFile = fopen(filename.c_str(), "wb");
    if (!vtkFile) {
        fmt::print(stderr, "Failed to open file {}\n", filename);
        return false;
    }
    fwrite(buffer.data(), 1, buffer.size(), vtkFile);
    fclose(vtkFile);

    return true;
}

// basic mesh interpolate
void Mesh3D::InterpolateTetraMesh(GRDECLData &data, int samples) {
    printf("Interpolate basic tetra mesh by hex mesh!\n");

    // 使用 OpenMP 并行区域，每个线程都有自己的随机数生成器
    #pragma omp parallel
    {
        // 每个线程创建一个自己的随机数生成器
        std::mt19937 gen(std::random_device{}());

        // 使用 OpenMP for 指令分发循环迭代
        #pragma omp for schedule(guided)
        for (INDEX i = 0; i < basic_tetra_size * 4; i += 4) { // only interpolate basic mesh
            INDEX corner_z, corner_y, corner_x;
            INDEX prev_y = 0, prev_x = 0;

            // 从 tetra mesh 中取出四个点
            Point3d p1 = points[tetrahedron_cells[i]];
            Point3d p2 = points[tetrahedron_cells[i + 1]];
            Point3d p3 = points[tetrahedron_cells[i + 2]];
            Point3d p4 = points[tetrahedron_cells[i + 3]];

            double tmp_young   = 0.0f;
            double tmp_poisson = 0.0f;
            double tmp_rockden = 0.0f;
            double tmp_poro    = 0.0f;

            // 对每个 tetrahedron 采样
            for (int s = 0; s < samples; s++) {
                // 采样 tetrahedron 内的随机点，使用当前线程的随机数生成器
                Point3d random_point = geo3d::samplePointInTetrahedron(p1, p2, p3, p4, gen);

                // 根据随机点在 hex mesh 中确定 cell 的索引
                data.GetCornerOfPoint(corner_z, corner_y, corner_x, random_point, prev_y, prev_x, gen);
                prev_y = corner_y;
                prev_x = corner_x;

                // 累加属性值
                tmp_young   += data.GetYoung(corner_z, corner_y, corner_x);
                tmp_poisson += data.GetPoisson(corner_z, corner_y, corner_x);
                tmp_rockden += data.GetRockDen(corner_z, corner_y, corner_x);
                tmp_poro    += data.GetPoro(corner_z, corner_y, corner_x);
            } // end samples

            // 计算平均值后写回 mesh
            young[i / 4]   = tmp_young   / samples;
            poisson[i / 4] = tmp_poisson / samples;
            rockden[i / 4] = tmp_rockden / samples;
            poro[i / 4]    = tmp_poro    / samples;
        } // end for
    } // end parallel

    printf("Interpolate basic tetra mesh done!\n");
}

// expand mesh interpolate
void Mesh3D::InterpolateExpandTetraMesh(GRDECLData &data, ExpandInfo &expand_info, int samples) {
    printf("Interpolate tetra mesh by expand mesh!\n");
    InterpolateTetraMesh(data, samples);

    // Process expand mesh
    if (!is_expand) {
        printf("No expand mesh, skip!\n");
        return;
    }

    bool positive_x, positive_y, positive_z;
    positive_x = data.GetCorner(0, 0, 0).coord[0][0][0].x < data.GetCorner(data.nz - 1, data.ny - 1, data.nx - 1).coord[1][1][1].x;
    positive_y = data.GetCorner(0, 0, 0).coord[0][0][0].y < data.GetCorner(data.nz - 1, data.ny - 1, data.nx - 1).coord[1][1][1].y;
    positive_z = data.GetCorner(0, 0, 0).coord[0][0][0].z < data.GetCorner(data.nz - 1, data.ny - 1, data.nx - 1).coord[1][1][1].z;

    // Get average value
    double average_young = data.avg_young;
    double average_poisson = data.avg_poisson;
    double average_rockden = data.avg_rockden;
    double average_poro = data.avg_poro;
    printf("average young: %f, poisson: %f, rockden: %f, poro: %f\n", average_young, average_poisson, average_rockden, average_poro);

    INDEX nz = data.nz;
    INDEX ny = data.ny;
    INDEX nx = data.nx;

    // preprocess: get 6 surfaces
    std::vector<Point2d> surface_z0;
    std::vector<Point2d> surface_z1;
    std::vector<Point2d> surface_y0;
    std::vector<Point2d> surface_y1;
    std::vector<Point2d> surface_x0;
    std::vector<Point2d> surface_x1;

    // convex judge
    std::vector<int> convex_z0;
    std::vector<int> convex_z1;
    std::vector<int> convex_y0;
    std::vector<int> convex_y1;
    std::vector<int> convex_x0;
    std::vector<int> convex_x1;

    // z0 surface
    surface_z0.reserve((ny + 1) * (nx + 1));
    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            surface_z0.push_back({data.GetCorner(0, j, i).coord[0][0][0].x, data.GetCorner(0, j, i).coord[0][0][0].y});
        }
        surface_z0.push_back({data.GetCorner(0, j, nx - 1).coord[0][0][1].x, data.GetCorner(0, j, nx - 1).coord[0][0][1].y});
    }
    for (INDEX i = 0; i < nx; i++) {
        surface_z0.push_back({data.GetCorner(0, ny - 1, i).coord[0][1][0].x, data.GetCorner(0, ny - 1, i).coord[0][1][0].y});
    }
    surface_z0.push_back({data.GetCorner(0, ny - 1, nx - 1).coord[0][1][1].x, data.GetCorner(0, ny - 1, nx - 1).coord[0][1][1].y});

    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            Point2d p1 = surface_z0[j * (nx + 1) + i];
            Point2d p2 = surface_z0[j * (nx + 1) + i + 1];
            Point2d p3 = surface_z0[(j + 1) * (nx + 1) + i];
            Point2d p4 = surface_z0[(j + 1) * (nx + 1) + i + 1];

            auto convex = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);
            if (convex.first) {
                convex_z0.push_back(0);
            } else {
                convex_z0.push_back(convex.second);
            }
        }
    }

    // z1 surface
    surface_z1.reserve((ny + 1) * (nx + 1));
    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            surface_z1.push_back({data.GetCorner(nz - 1, j, i).coord[1][0][0].x, data.GetCorner(nz - 1, j, i).coord[1][0][0].y});
        }
        surface_z1.push_back({data.GetCorner(nz - 1, j, nx - 1).coord[1][0][1].x, data.GetCorner(nz - 1, j, nx - 1).coord[1][0][1].y});
    }
    for (INDEX i = 0; i < nx; i++) {
        surface_z1.push_back({data.GetCorner(nz - 1, ny - 1, i).coord[1][1][0].x, data.GetCorner(nz - 1, ny - 1, i).coord[1][1][0].y});
    }
    surface_z1.push_back({data.GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].x, data.GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].y});

    for (INDEX j = 0; j < ny; j++) {
        for (INDEX i = 0; i < nx; i++) {
            Point2d p1 = surface_z1[j * (nx + 1) + i];
            Point2d p2 = surface_z1[j * (nx + 1) + i + 1];
            Point2d p3 = surface_z1[(j + 1) * (nx + 1) + i];
            Point2d p4 = surface_z1[(j + 1) * (nx + 1) + i + 1];

            auto convex = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);
            if (convex.first) {
                convex_z1.push_back(0);
            } else {
                convex_z1.push_back(convex.second);
            }
        }
    }

    // y0 surface
    surface_y0.reserve((nz + 1) * (nx + 1));
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX i = 0; i < nx; i++) {
            surface_y0.push_back({data.GetCorner(k, 0, i).coord[0][0][0].x, data.GetCorner(k, 0, i).coord[0][0][0].z});
        }
        surface_y0.push_back({data.GetCorner(k, 0, nx - 1).coord[0][0][1].x, data.GetCorner(k, 0, nx - 1).coord[0][0][1].z});
    }
    for (INDEX i = 0; i < nx; i++) {
        surface_y0.push_back({data.GetCorner(nz - 1, 0, i).coord[1][0][0].x, data.GetCorner(nz - 1, 0, i).coord[1][0][0].z});
    }
    surface_y0.push_back({data.GetCorner(nz - 1, 0, nx - 1).coord[1][0][1].x, data.GetCorner(nz - 1, 0, nx - 1).coord[1][0][1].z});

    for (INDEX k = 0; k < nz; k++) {
        for (INDEX i = 0; i < nx; i++) {
            Point2d p1 = surface_y0[k * (nx + 1) + i];
            Point2d p2 = surface_y0[k * (nx + 1) + i + 1];
            Point2d p3 = surface_y0[(k + 1) * (nx + 1) + i];
            Point2d p4 = surface_y0[(k + 1) * (nx + 1) + i + 1];

            auto convex = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);
            if (convex.first) {
                convex_y0.push_back(0);
            } else {
                convex_y0.push_back(convex.second);
            }
        }
    }

    // y1 surface
    surface_y1.reserve((nz + 1) * (nx + 1));
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX i = 0; i < nx; i++) {
            surface_y1.push_back({data.GetCorner(k, ny - 1, i).coord[0][1][0].x, data.GetCorner(k, ny - 1, i).coord[0][1][0].z});
        }
        surface_y1.push_back({data.GetCorner(k, ny - 1, nx - 1).coord[0][1][1].x, data.GetCorner(k, ny - 1, nx - 1).coord[0][1][1].z});
    }
    for (INDEX i = 0; i < nx; i++) {
        surface_y1.push_back({data.GetCorner(nz - 1, ny - 1, i).coord[1][1][0].x, data.GetCorner(nz - 1, ny - 1, i).coord[1][1][0].z});
    }
    surface_y1.push_back({data.GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].x, data.GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].z});

    for (INDEX k = 0; k < nz; k++) {
        for (INDEX i = 0; i < nx; i++) {
            Point2d p1 = surface_y1[k * (nx + 1) + i];
            Point2d p2 = surface_y1[k * (nx + 1) + i + 1];
            Point2d p3 = surface_y1[(k + 1) * (nx + 1) + i];
            Point2d p4 = surface_y1[(k + 1) * (nx + 1) + i + 1];

            auto convex = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);
            if (convex.first) {
                convex_y1.push_back(0);
            } else {
                convex_y1.push_back(convex.second);
            }
        }
    }

    // x0 surface
    surface_x0.reserve((nz + 1) * (ny + 1));
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            surface_x0.push_back({data.GetCorner(k, j, 0).coord[0][0][0].y, data.GetCorner(k, j, 0).coord[0][0][0].z});
        }
        surface_x0.push_back({data.GetCorner(k, ny - 1, 0).coord[0][1][0].y, data.GetCorner(k, ny - 1, 0).coord[0][1][0].z});
    }
    for (INDEX j = 0; j < ny; j++) {
        surface_x0.push_back({data.GetCorner(nz - 1, j, 0).coord[1][0][0].y, data.GetCorner(nz - 1, j, 0).coord[1][0][0].z});
    }
    surface_x0.push_back({data.GetCorner(nz - 1, ny - 1, 0).coord[1][1][0].y, data.GetCorner(nz - 1, ny - 1, 0).coord[1][1][0].z});

    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            Point2d p1 = surface_x0[k * (ny + 1) + j];
            Point2d p2 = surface_x0[k * (ny + 1) + j + 1];
            Point2d p3 = surface_x0[(k + 1) * (ny + 1) + j];
            Point2d p4 = surface_x0[(k + 1) * (ny + 1) + j + 1];

            auto convex = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);
            if (convex.first) {
                convex_x0.push_back(0);
            } else {
                convex_x0.push_back(convex.second);
            }
        }
    }

    // x1 surface
    surface_x1.reserve((nz + 1) * (ny + 1));
    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            surface_x1.push_back({data.GetCorner(k, j, nx - 1).coord[0][0][1].y, data.GetCorner(k, j, nx - 1).coord[0][0][1].z});
        }
        surface_x1.push_back({data.GetCorner(k, ny - 1, nx - 1).coord[0][1][1].y, data.GetCorner(k, ny - 1, nx - 1).coord[0][1][1].z});
    }
    for (INDEX j = 0; j < ny; j++) {
        surface_x1.push_back({data.GetCorner(nz - 1, j, nx - 1).coord[1][0][1].y, data.GetCorner(nz - 1, j, nx - 1).coord[1][0][1].z});
    }
    surface_x1.push_back({data.GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].y, data.GetCorner(nz - 1, ny - 1, nx - 1).coord[1][1][1].z});

    for (INDEX k = 0; k < nz; k++) {
        for (INDEX j = 0; j < ny; j++) {
            Point2d p1 = surface_x1[k * (ny + 1) + j];
            Point2d p2 = surface_x1[k * (ny + 1) + j + 1];
            Point2d p3 = surface_x1[(k + 1) * (ny + 1) + j];
            Point2d p4 = surface_x1[(k + 1) * (ny + 1) + j + 1];

            auto convex = geo2d::checkQuadrilateralConvexity(p1, p2, p3, p4);
            if (convex.first) {
                convex_x1.push_back(0);
            } else {
                convex_x1.push_back(convex.second);
            }
        }
    }



    // Set reservoir cell data
    for (INDEX i = 0; i < basic_tetra_size; i++) {
        cell_types[i] = ExpandCellType::RESERVOIR;
    }

    // traverse expand mesh
    #pragma omp parallel for schedule(guided)
    for (INDEX i = basic_tetra_size * 4; i < expand_tetra_size * 4; i += 4) {
        INDEX corner_z = -1;
        INDEX corner_y = -1;
        INDEX corner_x = -1;
        INDEX inter_z = -1;
        INDEX inter_y = -1;
        INDEX inter_x = -1;

        Point3d p1 = points[tetrahedron_cells[i]];
        Point3d p2 = points[tetrahedron_cells[i + 1]];
        Point3d p3 = points[tetrahedron_cells[i + 2]];
        Point3d p4 = points[tetrahedron_cells[i + 3]];

        Point3d center_pos = expand_info.GetTetrahedronPosition(p1, p2, p3, p4); // 3x3x3 value
        Point3d center_point = geo3d::computeMasscenter(p1, p2, p3, p4);

        if ((positive_z && center_pos.z == 1) || (!positive_z && center_pos.z == -1)) {
            corner_z = data.nz - 1;
        }
        else if ((positive_z && center_pos.z == -1) || (!positive_z && center_pos.z == 1)) {
            corner_z = 0;
        }
        
        if ((positive_y && center_pos.y == 1) || (!positive_y && center_pos.y == -1)) {
            corner_y = data.ny - 1;
        }
        else if ((positive_y && center_pos.y == -1) || (!positive_y && center_pos.y == 1)) {
            corner_y = 0;
        }

        if ((positive_x && center_pos.x == 1) || (!positive_x && center_pos.x == -1)) {
            corner_x = data.nx - 1;
        }
        else if ((positive_x && center_pos.x == -1) || (!positive_x && center_pos.x == 1)) {
            corner_x = 0;
        }
        
        // cant interpolate poro, directly use average value
        poro[i / 4] = average_poro;

        // TODO: use function
        if (corner_z == -1) { // middle layer
            cell_types[i / 4] = ExpandCellType::SIDESBURDEN;
            if (corner_y != -1 && corner_x != -1) { // z axis search
                inter_x = corner_x;
                inter_y = corner_y;
                int xval = (corner_x == 0) ? 0 : 1;
                int yval = (corner_y == 0) ? 0 : 1;
                for (int i = 0; i < data.nz; i++) {
                    double z0 = data.GetCorner(i, corner_y, corner_x).coord[0][yval][xval].z;
                    double z1 = data.GetCorner(i, corner_y, corner_x).coord[1][yval][xval].z; 
                    if ((center_point.z - z0) * (center_point.z - z1) <= 0) {
                        inter_z = i;
                        break;
                    }
                }
            }
            else if (corner_y != -1) { // y plane search
                inter_y = corner_y;
                if (corner_y == 0) {
                    geo2d::searchPointOnSurface(surface_y0, convex_y0, {center_point.x, center_point.z}, inter_x, inter_z, nx, nz);
                }
                else if (corner_y == ny - 1) {
                    geo2d::searchPointOnSurface(surface_y1, convex_y1, {center_point.x, center_point.z}, inter_x, inter_z, nx, nz);
                }
            }
            else if (corner_x != -1) { // x plane search
                inter_x = corner_x;
                if (corner_x == 0) {
                    geo2d::searchPointOnSurface(surface_x0, convex_x0, {center_point.y, center_point.z}, inter_y, inter_z, ny, nz);
                }
                else if (corner_x == nx - 1) {
                    geo2d::searchPointOnSurface(surface_x1, convex_x1, {center_point.y, center_point.z}, inter_y, inter_z, ny, nz);
                }
            }
            else { // basic mesh, invalid
                printf("Error: invalid expand cell\n");
            }
        }
        else if (corner_z == 0) { // bottom layer
            inter_z = corner_z;
            if (corner_y == -1 && corner_x == -1) { // z bottom plane search
                cell_types[i / 4] = ExpandCellType::UNDERBURDEN;
                geo2d::searchPointOnSurface(surface_z0, convex_z0, {center_point.x, center_point.y}, inter_x, inter_y, nx, ny);
            }
            else {
                cell_types[i / 4] = ExpandCellType::UNDERSIDES;
                if (corner_y != -1 && corner_x != -1) { // 角点直接赋值
                    inter_x = corner_x;
                    inter_y = corner_y;
                }
                else if (corner_y != -1) { // x axis search
                    inter_y = corner_y;
                    int yval = (corner_y == 0) ? 0 : 1;
                    int zval = 0;
                    for (int i = 0; i < data.nx; i++) {
                        double x0 = data.GetCorner(corner_z, corner_y, i).coord[zval][yval][0].x;
                        double x1 = data.GetCorner(corner_z, corner_y, i).coord[zval][yval][1].x;
                        if ((center_point.x - x0) * (center_point.x - x1) <= 0) {
                            inter_x = i;
                            break;
                        }
                    }
                }
                else if (corner_x != -1) { // y axis search
                    inter_x = corner_x;
                    int xval = (corner_x == 0) ? 0 : 1;
                    int zval = 0;
                    for (int i = 0; i < data.ny; i++) {
                        double y0 = data.GetCorner(corner_z, i, corner_x).coord[zval][0][xval].y;
                        double y1 = data.GetCorner(corner_z, i, corner_x).coord[zval][1][xval].y;
                        if ((center_point.y - y0) * (center_point.y - y1) <= 0) {
                            inter_y = i;
                            break;
                        }
                    }
                }
            }
        }
        else { // corner_z == data.nz - 1, top layer
            inter_z = corner_z;
            if (corner_y == -1 && corner_x == -1) { // z top plane search
                cell_types[i / 4] = ExpandCellType::OVERBURDEN;
                geo2d::searchPointOnSurface(surface_z1, convex_z1, {center_point.x, center_point.y}, inter_x, inter_y, nx, ny);
            }
            else {
                cell_types[i / 4] = ExpandCellType::OVERSIDES;
                if (corner_y != -1 && corner_x != -1) { // 角点直接赋值
                    inter_x = corner_x;
                    inter_y = corner_y;
                }
                else if (corner_y != -1) { // x axis search
                    inter_y = corner_y;
                    int yval = (corner_y == 0) ? 0 : 1;
                    int zval = 1;
                    for (int i = 0; i < data.nx; i++) {
                        double x0 = data.GetCorner(corner_z, corner_y, i).coord[zval][yval][0].x;
                        double x1 = data.GetCorner(corner_z, corner_y, i).coord[zval][yval][1].x;
                        if ((center_point.x - x0) * (center_point.x - x1) <= 0) {
                            inter_x = i;
                            break;
                        }
                    }
                }
                else if (corner_x != -1) { // y axis search
                    inter_x = corner_x;
                    int xval = (corner_x == 0) ? 0 : 1;
                    int zval = 1;
                    for (int i = 0; i < data.ny; i++) {
                        double y0 = data.GetCorner(corner_z, i, corner_x).coord[zval][0][xval].y;
                        double y1 = data.GetCorner(corner_z, i, corner_x).coord[zval][1][xval].y;
                        if ((center_point.y - y0) * (center_point.y - y1) <= 0) {
                            inter_y = i;
                            break;
                        }
                    }
                }
            }
        }

        // interpolate get young, poisson, rockden
        // if (inter_z == -1 || inter_y == -1 || inter_x == -1) {
            // printf("Error: inter_z, inter_y, inter_x not found\n");
        // }

        young[i / 4] = data.GetYoung(inter_z, inter_y, inter_x);
        poisson[i / 4] = data.GetPoisson(inter_z, inter_y, inter_x);
        rockden[i / 4] = data.GetRockDen(inter_z, inter_y, inter_x);
    }
    printf("Interpolate expand tetra mesh done!\n");
}
